Szymon Malec, Viktoriia Morhun
Ruch uliczny wzrasta w ostatnich latach bardzo szybko, co potwierdzają wyniki Generalnego Pomiaru Ruchu przeprowadzanego co pięć lat przez GDDKiA. W ciągu ostatnich pięciu lat średni dobowy ruch w Polsce wzrósł aż o 21% [1]. Dynamicznie rozwijające się miasta takie, jak Wrocław, gdzie cały czas powstają nowe osiedla i przybywa mieszkańców, muszą radzić sobie ze wzrastającym natężeniem na drogach. To zmusza władze do budowania nowych dróg, rozwijania komunikacji miejskiej (by odciągnąć ludzi od samochodów) czy wprowadzania inteligentnych systemów sygnalizacji świetlnej. Postanowiliśmy przyjrzeć się problemowi oraz prawom, którymi rządzą się korki drogowe, by następnie przenieść nasze obserwacje na równania matematyczne w celu przewidzenia długości takich korków. Tak udało nam się sformułować dwa modele, które następnie zastosowaliśmy do jednej z wrocławskich ulic.
Zanim wyprowadzimy równanie opisujące długość korka, zastanówmy się, kiedy powstaje korek drogowy i jakie są tego przyczyny. Potrzebne będzie nam do tego kilka podstawowych pojęć związanych z tematem:
natężenie ruchu ($n$) - iloczyn ilości aut $a$ przejężdżających przez dany punkt na drodze o swobodnym ruchu w jednostce czasu i średniej długości samochodu (włącznie z bezpiecznym odstępem) $r$. Natężenie możemy otrzymać licząc ilość przejeżdżających aut $a$ przez pewien czas $t$, a następnie podstawiąjąc do wzoru $$ n = \frac{ar}{t}, $$
gęstość ruchu ($g$) - część drogi jaka jest zajęta przez jadące samochody na danym odcinku w jednostce czasu. Gęstość możemy otrzymać po zmierzeniu natężenia $n$, prędkości średniej $v_{śr}$ jadących aut oraz średniej długości samochodu (włącznie z bezpiecznym odstępem) $r$, a następnie korzystając ze wzoru $$ g = \frac{n}{v_{śr}} = \frac{ar}{v_{śr}t}, $$ przy czym zakładamy, że $g \in [0, 1)$, ponieważ gdyby $g$ było równe 1, to przez bardzo małe odstępy auta nie byłyby w stanie się poruszać,
przepustowość ($p$) - iloczyn ilości aut $a$ omijających punkt spowalniający ruch na drodze oraz średniej długości samochodu (włącznie z bezpiecznym odstępem) $r$ w jednostce czasu. Żeby otrzymać przepustowość, korzystamy ze wzoru $$ p = \frac{ar}{t}. $$
Korek tworzy się wtedy, gdy w jakimś punkcie na drodze znajduje się lub pojawia się coś, co spowalnia ruch (zmniejsza przepustowość w tym punkcie), przy czym przepustowość musi być mniejsza od natężenia tzn. ilość aut dojeżdżających musi być większa od aut odjeżdżających.

Przykłady miejsc/sytuacji, w których tworzą się korki:
Do ułożenia równania opisującego długość korka w czasie przydatne okazuje się być porównanie samochodów do płynu. Koncept gęstości ruchu i porównania korka do płynu zaczerpnęliśmy z [2]. By było to możliwe, musimy założyć, że wszystkie auta w korku ruszają jednocześnie, bez opóźnienia. Wtedy korek będzie zachowywał się podobnie do płynu. Drogę interpretujemy jako pionowy kanał, którym płynie woda z prędkością $v_{śr}$. W pewnym momencie kanał ten zostaje częściowo zatkany. Gęstość w tym przypadku oznacza procent wypełnienia kanału wodą podczas swobodnego przepływu, natężenie - ilość dopływającej wody, a przepustowość - ilość wody, która przepływa przez miejsce zatkania. Długością korka będzie poziom wody w kanale od miejsca zatkania - oznaczmy $L$. Jeśli $L(t)$ będzie opisywało poziom wody (długość korka) w czasie, to $L'(t)$ będzie prędkością z jaką ten poziom (korek) rośnie.

Ponieważ prędkość płynącej wody $v_{śr}$ oraz prędkość poziomu wody $L'$ mają przeciwne zwroty, to przyrost tego poziomu spowodowany dopłynięciem nowej wody w momencie $t$ będzie wynosił $ \left( v_{śr} + L' \right) g(t) $. Odejmując od tego przepustowość $p(t)$, otrzymamy równanie na przyrost poziomu wody: $$ L' = \left( v_{śr} + L' \right) g(t) - p(t). $$ Tak przy pomocy odpowiedniego porównania udało nam się otrzymać równanie różniczkowe opisujące długość korka. Następnie po przekształceniu otrzymamy $$ L' = \frac{v_{śr} g(t) - p(t)}{1 - g(t)}. $$ Rozwiązując równanie dostajemy $$ L(t) = L_{0} + \int_{0}^{t} \frac{v_{śr} g(t) - p(t)}{1 - g(t)} \mathrm{d}t, $$ gdzie $L_{0}$ - początkowy korek.
Istotnym założeniem jest tutaj to, że $L \geq 0$, co oznacza, że nie każde równanie będzie dało się rozwiązać analitycznie. Będzie to możliwe tylko w przypadku, gdy korek utrzymuje się cały czas i jego długość nie spada poniżej zera. W innym przypadku trzeba wykorzystać metody numeryczne.
Poprzednie równanie opisywało długość korka przy założeniu, że wszystkie auta ruszają jednocześnie. Z doświadczenia wiemy, że w rzeczywistości jest inaczej. Auta ruszają jedno za drugim z pewnym opóźnieniem, które nazywać będziemy czasem reakcji. W wyniku tego powstaje zjawisko przypominające falę, którą z kolei nazwiemy falą drogową.

Żeby lepiej zrozumieć to zjawisko rozpatrzmy przykład. Z powodu czerwonych świateł utworzył się korek. Gdy pojawia się zielone światło, pierwsze auto rusza, następnie po czasie reakcji rusza kolejne itd. Powstaje wcześniej wspomniana fala, która porusza się z jakąś prędkością $v_f$. Ostatnie auto ruszy dopiero, gdy fala ta dotrze do niego, czyli pokona odległość równą długości korka. Widzimy zatem, że korek skróci się po pewnym czasie $\Delta t$ od zapalenia się zielonych świateł.
Ponieważ czas $\Delta t$ to czas w jakim fala drogowa poruszająca się ze stałą predkością $v_f$ pokona długość korka $L$, to $$\Delta t = \frac{L}{v_f}.$$ Wtedy przepustowość w momencie $t$ będzie równa $p(t - \Delta t) = p\left(t - \frac{L}{v_f}\right)$.
Rozważmy teraz przypadek, w którym w momencie $t = 0$ utworzył się korek o długości $L_0$ oraz natężenie od tego momentu jest zerowe. Mamy jedynie przepustowość $p(t - \frac{L}{v_f})$. Wtedy długość korka można opisać równaniem $$ L(t) = L_0 - \int_{0}^{\gamma}p(\tau) \mathrm{d}\tau, $$ gdzie $\gamma = t - \frac{L}{v_f}$. Po zróżniczkowaniu obustronnie względem $\gamma$ otrzymujemy $$ \frac{\mathrm{d}L}{\mathrm{d}\gamma} = - p(\gamma). $$ Następnie domnażamy $\frac{\mathrm{d}\gamma}{\mathrm{d}t}$ i upraszczamy: $$ \frac{\mathrm{d}L}{\mathrm{d}\gamma} \frac{\mathrm{d}\gamma}{\mathrm{d}t} = - p(\gamma) \frac{\mathrm{d}\gamma}{\mathrm{d}t}, $$ $$ \frac{\mathrm{d}L}{\mathrm{d}t} = - p(\gamma) \frac{\mathrm{d}}{\mathrm{d}t} \left(t - \frac{L}{v_f}\right), $$ $$ L' = - \left(1 - \frac{L'}{v_f}\right) p\left(t - \frac{L}{v_f}\right). $$ W ten spośob otrzymaliśmy postać, jaką ma przepustowość w modelu z opóźnieniem. Teraz możemy sformułować pełne równanie zawierające natężenie oraz przepustowość. Ponieważ opóźnienie $\Delta t$ nie wpływa na przyrost nowych samochodów, będzie on taki sam, jak w pierwszym modelu tj. $ \left( v_{śr} + L' \right) g(t) $. Zatem równanie opisujące długość korka będzie mieć postać $$ L' = \left( v_{śr} + L' \right) g(t) - \left(1 - \frac{L'}{v_f}\right) p\left(t - \frac{L}{v_f} \right). $$ Po przekształceniu otrzymujemy $$ L' = \frac{v_{śr} g(t) - p\left(t - \frac{L}{v_f}\right)}{1 - g(t) - \frac{1}{v_f} p \left(t - \frac{L}{v_f} \right)}. $$ Należy tutaj zaznaczyć, że poza założeniem $L \geq 0$, dochodzi tutaj jeszcze jedno istotne założenie tj. $L' \leq v_f$ wynikające stąd, że gdyby $L' > v_f$, to $$ \left( 1 - \frac{L'}{v_f} \right) < 0, $$ a wtedy $$ \left(1 - \frac{L'}{v_f}\right) p\left(t - \frac{L}{v_f} \right) < 0, $$ co byłoby sprzeczne, ponieważ przepustowość nie może być ujemna.
Najważniejszą częścią naszego projektu było zastosowanie przedstawionych modeli do prawdziwej sytuacji. Przyjrzeliśmy się zatem najbardziej zakorkowanym ulicom we Wrocławiu. Miejscem, które wydało nam się najbardziej odpowiednie do zamodelowania jest most Pokoju oraz fragment ulicy Wyszyńskiego. Na poniższych zdjęciach możemy zobaczyć, jak wygląda wspomniane miejsce oraz który odcinek badamy.

Długość badanej drogi wynosi ok. 416m, co zmierzyliśmy za pomocą Map Google. Droga składa się z dwóch głównych pasów oraz z dwóch dodatkowych (przeznaczonych do jazdy w prawo i w lewo), które znajdują się już przy samych światłach. Na całej długości odcinka nie pojawia się żaden zjazd ani nic, co mogłoby spowalniać ruch (np. przejście dla pieszych). Czas, przez który świeci się światło zielone i czerwone jest nieregularny ze względu na zastosowany we Wrocławiu system ITS (Intelligent Transport System), który analizuje ruch uliczny i dopasowuje czas świateł do sytuacji. Postanowiliśmy zatem zbadać, jaki stosunek musi zachować światło zielone do czerwonego, by długość korka nie przekroczyła długości badanego odcinka. Aby nasz model jak najlepiej oddawał rzeczywistość, wykonaliśmy serię niezbędnych pomiarów.
# Wczytujemy potrzebne biblioteki
using Plots
using LsqFit
using DataTables
By otrzymać prędkość średnią, mierzyliśmy czas jaki zajmowało samochodowi pokonanie odcinka od pasów na początku drogi do początku mostu. Następnie podzieliliśmy długość tego odcinka, która wynosi 114m, przez średni czas i dostaliśmy $v_{śr} \approx 11,12\frac{\mathrm{m}}{\mathrm{s}}$.
times = DataTable(Czas = [8.42, 6.77, 10.09, 7.12, 7.30, 7.88, 9.02, 7.35, 7.28, 6.51, 9.87, 7.41, 7.57, 9.65, 8.69,
5.80, 6.52, 7.49, 9.91, 10.17, 7.11, 7.56, 5.78, 5.77, 6.84, 7.23, 7.86, 9.17, 8.73, 10.52])
average_time = sum(k[1] for k in times[:, 1]) / length(times[:, 1])
v_śr = 88 / average_time
println("Prędkość średnia: ", v_śr)
Prędkość średnia: 11.12094022494629
W celu odnalezienia długości średniej samochodu wraz z odstępem, liczyliśmy ile aut stojących w korku na jednym z pasów mieści się na odcinku od świateł do mostu, po czym podzieliliśmy długość tego odcinka (180m) przez średnią liczbę aut.
cars = DataTable(Auta = [25, 23, 26, 27, 25, 26, 24, 25])
average_cars = sum(k[1] for k in cars[:, 1]) / length(cars[:, 1])
r_śr = 180 / average_cars
r = r_śr/2
show(cars)
println("\n\nŚrednia długość samochodu: ", r_śr)
println("r = ", r)
8x1 DataTable Auta ──── 25 23 26 27 25 26 24 25 Średnia długość samochodu: 7.164179104477612 r = 3.582089552238806
Otrzymany wynik dzielimy jeszcze przez dwa ze względu na dwa pasy ruchu, stąd $r \approx 3,58\mathrm{m}$.
Aby znaleźć prędkość tej fali, mierzyliśmy czas od momentu ruszenia pierszego samochodu w korku do momentu ruszenia ostatniego samochodu znajdującego się na odcinku od świateł do mostu. Następnie z otrzymanych pomiarów wyliczyliśmy czas średni i podzieliliśmy długość odcinka przez ten czas.
times = DataTable(Czas = [26.30, 30.39, 27.54, 30.78, 27.73, 31.44, 31.57, 26.23, 24.59, 29.14, 28.58])
average_time = sum(k[1] for k in times[:, 1]) / length(times[:, 1])
v_f = 180 / average_time
show(times)
println("\n\nPrędkość fali: ", v_f)
11x1 DataTable Czas ───── 26.3 30.39 27.54 30.78 27.73 31.44 31.57 26.23 24.59 29.14 ... with 1 more row Prędkość fali: 6.299914092080565
A zatem $v_f \approx 6,3 \frac{\mathrm{m}}{\mathrm{s}}$.
Gęstość zmierzyliśmy licząc ilość aut wjeżdżających na badany odcinek w ciągu kolejnych 20-minututowych okresów czasu od godz. 14.00 do 18.00. Wybraliśmy te godziny, ponieważ jest to pora, o której ma miejsce największy ruch, a miasto jest najbardziej zakorkowane. Następnie, korzystając ze wzoru na gęstość, otrzymaliśmy poniższe dane:
T = 14400
a_data = [443, 459, 460, 472, 502, 488, 493, 464, 442, 418, 421, 381]
g_data = (a_data .* r) ./ (1200 * v_śr)
scatter(1/3600+14:1200/3600:T/3600+14, g_data, xlabel="Godzina", ylabel="g", legend=false)
Korzystając z biblioteki LsqFit dostępnej w Julii, udało nam się dopasować krzywą do danych. Jako funkcję bazową przyjęliśmy $$ g(t) = \alpha \cdot \mathrm{sin}(\omega t + \beta) + \delta, $$ gdzie $\alpha, \omega, \beta, \delta$ są szukanymi parametrami.
g_mod(t, p) = p[1] .* sin.(p[2] .* t .+ p[3]) .+ p[4]
g_fit = curve_fit(g_mod, 1:1200:T, g_data, [1.0, 0.0005, 0.0, 0.0])
scatter(1/3600+14:1200/3600:T/3600+14, g_data, xlabel="Godzina", ylabel="g", label=false)
plt = plot!(1/3600+14:1/3600:T/3600+14, g_mod(1:T, g_fit.param), label="g(t)")
display(plt)
symbols = ["α", "ω", "β", "δ"]
for param in zip(symbols, g_fit.param)
println(param[1], " = ", param[2])
end
α = 0.017363738729266526 ω = 0.00027397276238646135 β = 0.15374465556085215 δ = 0.11422808354362937
Z wykresu widzimy jak gęstość zmienia się w godzinach szczytu.
Przepustowość otrzymaliśmy licząc ilość aut przejeżdżających na zielonym świetle, dzieląc przez ten czas oraz domnażając $r$ (według wzoru na $p$).
capacity = DataTable(Czas = [32, 36, 34, 47, 47, 38, 32, 44, 31, 42, 49, 37, 36, 38], Auta = [35, 39, 35, 51, 50, 36, 32, 46, 30, 46, 50, 40, 37, 40])
average_capacity = r * sum([k[2]/k[1] for k in capacity[:,1]]) / length(capacity[:,1])
show(capacity)
print("\n\nŚrednia przepustowość: ", average_capacity)
14x2 DataTable Czas │ Auta ─────┼───── 32 │ 35 36 │ 39 34 │ 35 47 │ 51 47 │ 50 38 │ 36 32 │ 32 44 │ 46 31 │ 30 42 │ 46 ... with 4 more rows Średnia przepustowość: 3.7338506268365004
Otrzymany wynik dotyczy oczywiście przepustowości podczas światła zielonego. Podczas gdy świeci się światło czerwone, wynosi ona zero.
Załóżmy, że łączny czas światła zielonego i czerwonego to 100s. Niech $T_Z$ - czas światła zielonego, $T_C$ - czas światła czerwonego, czyli $T_Z + T_C = 100\mathrm{s}$.
Wyniki pomiarów podstawiliśmy do pierwszego modelu, czyli
$$ L' = \frac{v_{śr} g(t) - p(t)}{1 - g(t)}, $$
gdzie
$$ g(t) = \alpha \cdot \mathrm{sin}(\omega t + \beta) + \delta $$
($\alpha, \omega, \beta, \delta$ takie jak dostaliśmy w 4.2.4),
$$ p(t) = \left\{ \begin{array}{lrl} 3,734 & \mathrm{dla} & 100n < t \leq 100n + T_Z \\ 0 & \mathrm{dla} & 100n + T_Z < t \leq 100(n+1) \\ \end{array} \right. , \hspace{0.5cm} n \in \mathbb{N}. $$
Spodziewamy się, że długość korka może spaść poniżej zera, a więc, żeby uniknąć sprzeczności, do wyliczenia rozwiązania skorzystamy z metod numerycznych, a konkretnie z metody Eulera niejawnej, czyli
$$ L_n = L_{n-1} + a \cdot L'(t_{n}), $$
gdzie
$t_n = 0,1n [ \mathrm{s} ],\ \ \ n = 1, 2, ...$
$a = t_n - t_{n-1} = 0,1\mathrm{s}$.
Pryjmijmy $L_0 = 0$. W kodzie uwzględnimy warunek, że $L \geq 0$.
Rozwiążemy równanie dla różnych wartości $T_Z$, by sprawdzić jaki musi być minimalny czas światła zielonego, by długość korka nie przekroczyła długości drogi (416m).
g(t) = g_fit.param[1] * sin(g_fit.param[2] * t + g_fit.param[3]) + g_fit.param[4]
function p(t, T_Z)
if 0 < t % 100 <= T_Z
return 3.734
else
return 0.0
end
end
a = 0.1 # Krok czasowy
anim1 = @animate for T_Z in 36:0.1:40
# Wartości L'(t)
dL = a * (v_śr .* g.(a:a:T) .- p.(a:a:T, T_Z)) ./ (1 .- g.(a:a:T))
# Metodą Eulera wyliczamy wartości L(t)
L = [dL[1]]
for t in 2a:a:T
index = Int(round(t/a))
l = L[index-1] + dL[index]
if l <= 0
append!(L, 0)
else
append!(L, l)
end
end
plot(a/3600+14:a/3600:T/3600+14, L, ylims=(0, 1000), title="T_z = $T_Z", xlabel="Godzina", ylabel="L", label="L(t)")
plot!(x->x, x->416, a/3600+14:a/3600:T/3600+14, label="416")
end
gif(anim1, "Korek1.gif", fps= 10)
┌ Info: Saved animation to │ fn = C:\Users\szymo\korki\Korek1.gif └ @ Plots C:\Users\szymo\.julia\packages\Plots\FKcum\src\animation.jl:104
Sprawdźmy, jak wyglądają wykresy dla $T_Z = 37.4\mathrm{s}$ i $T_Z = 37.5\mathrm{s}$.
T_Z1 = 37.4
dL = a * (v_śr .* g.(a:a:T) .- p.(a:a:T, T_Z1)) ./ (1 .- g.(a:a:T))
L = [dL[1]]
for t in 2a:a:T
index = Int(round(t/a))
l = L[index-1] + dL[index]
if l <= 0
append!(L, 0)
else
append!(L, l)
end
end
plot(a/3600+14:a/3600:T/3600+14, L, title="T_z = $T_Z1", xlabel="Godzina", ylabel="L", label="L(t)")
plt1 = plot!(x->x, x->416, a/3600+14:a/3600:T/3600+14, label="416")
T_Z2 = 37.5
dL = a * (v_śr .* g.(a:a:T) .- p.(a:a:T, T_Z2)) ./ (1 .- g.(a:a:T))
L = [dL[1]]
for t in 2a:a:T
index = Int(round(t/a))
l = L[index-1] + dL[index]
if l <= 0
append!(L, 0)
else
append!(L, l)
end
end
plot(a/3600+14:a/3600:T/3600+14, L, title="T_z = $T_Z2", xlabel="Godzina", ylabel="L", label="L(t)")
plt2 = plot!(x->x, x->416, a/3600+14:a/3600:T/3600+14)
plot(plt1, plt2, layout=(1, 2), size=(1000, 380))
Widzimy, że dla $T_Z = 37.4$ korek wykroczył już poza ulicę. Zatem według tego modelu światło zielone musi świecić się przez przynajmniej 37,5s. Czyli stosunek światła zielonego do czerwonego wynosi $$ \frac{T_Z}{T_C} = \frac{37,5}{62,5} = 0,6. $$
Przypomnijmy równanie tego modelu: $$ L' = \frac{v_{śr} g(t) - p\left(t - \frac{L}{v_f}\right)}{1 - g(t) - \frac{1}{v_f} p \left(t - \frac{L}{v_f} \right)}. $$ Funkcje $g$ i $p$ są takie, jak w poprzednim podpunkcie.
Do rozwiązania tego równania skorzystaliśmy z metody Eulera jawnej, czyli
$$ L_n = L_{n-1} + a \cdot \frac{v_{śr} \cdot g(t_{n-1}) - p\left(t_{n-1} - \frac{L_{n-1}}{v_f}\right)}{1 - g(t_{n-1}) - \frac{1}{v_f} p \left(t_{n-1} - \frac{L_{n-1}}{v_f} \right)}, $$
gdzie
$t_n = 0,1n [ \mathrm{s} ],\ \ \ n = 1, 2, ...$
$a = t_n - t_{n-1} = 0,1\mathrm{s}$.
Ponownie przyjmijmy $L_0 = 0$ i rozwiążmy równanie dla różnych wartości $T_Z$.
a = 0.1 # Krok czasu
anim2 = @animate for T_Z in 36:0.1:40
L = [0.0]
for t in a:a:T-a
index = Int(round(t/a))
dL = a * (v_śr * g(t) - p(t - L[index]/v_f, T_Z)) / (1 - g(t) - p(t - L[index]/v_f, T_Z)/v_f)
l = L[index] + dL
if l <= 0
append!(L, 0)
else
append!(L, l)
end
end
plot(a/3600+14:a/3600:T/3600+14, L, ylims=(0, 1100), title="T_z = $T_Z", xlabel="Godzina", ylabel="L", label="L(t)")
plot!(x -> x, x -> 416, a/3600+14:a/3600:T/3600+14, label="416")
end
gif(anim2, "Korek2.gif", fps= 10)
┌ Info: Saved animation to │ fn = C:\Users\szymo\korki\Korek2.gif └ @ Plots C:\Users\szymo\.julia\packages\Plots\FKcum\src\animation.jl:104
Sprawdźmy, jak wyglądają wykresy dla $T_Z = 37.7\mathrm{s}$ i $T_Z = 37.8\mathrm{s}$.
T_Z1 = 37.7
L = [0.0]
for t in a:a:T-a
index = Int(round(t/a))
dL = a * (v_śr * g(t) - p(t - L[index]/v_f, T_Z1)) / (1 - g(t) - p(t - L[index]/v_f, T_Z1)/v_f)
l = L[index] + dL
if l <= 0
append!(L, 0)
else
append!(L, l)
end
end
plot(a/3600+14:a/3600:T/3600+14, L, title="T_z = $T_Z1", xlabel="Godzina", ylabel="L[m]", label="L(t)")
plt1 = plot!(x -> x, x -> 416, a/3600+14:a/3600:T/3600+14, label="416")
T_Z2 = 37.8
L = [0.0]
for t in a:a:T-a
index = Int(round(t/a))
dL = a * (v_śr * g(t) - p(t - L[index]/v_f, T_Z2)) / (1 - g(t) - p(t - L[index]/v_f, T_Z2)/v_f)
l = L[index] + dL
if l <= 0
append!(L, 0)
else
append!(L, l)
end
end
plot(a/3600+14:a/3600:T/3600+14, L, title="T_z = $T_Z2", xlabel="Godzina", ylabel="L", label="L(t)")
plt2 = plot!(x -> x, x -> 416, a/3600+14:a/3600:T/3600+14, label="416")
plot(plt1, plt2, layout=(1, 2), size=(1000, 380))
Widzimy, że dla $T_Z = 37.7$ korek wykroczył już poza ulicę. Zatem według modelu drugiego światło zielone musi świecić się przez przynajmniej 37,8s. Czyli stosunek światła zielonego do czerwonego wynosi $$ \frac{T_Z}{T_C} = \frac{37,8}{62,2} \approx 0,6077. $$
Porównując obydwa modele, zauważamy, że ich wyniki różnią się nieznacznie, co oznacza, że fala drogowa nie ma w tym przypadku znaczącego wpływu na długość korka. W obu przypadkach stosunek światła zielonego do czerwonego wynosi w przybliżeniu 3:5. Żeby porównać nasze wyniki z rzeczywistością, wykonaliśmy kilka pomiarów czasu świateł.
lights = DataTable(Czerwone = [60, 57, 54, 62, 65, 59, 69, 51], Zielone = [32, 36, 34, 47, 47, 32, 44, 31])
show(lights)
green = sum([k[2] for k in lights[:,1]]) / length(lights[:,1])
red = sum([k[1] for k in lights[:,1]]) / length(lights[:,1])
green_red_ratio = green/red
println("\n\nŚredni czas światła zielonego: ", green)
println("Średni czas światła czerwonego: ", red)
println("Stosunek czasu swiatła zielonego do czerwonego: ", green_red_ratio)
8x2 DataTable Czerwone │ Zielone ─────────┼──────── 60 │ 32 57 │ 36 54 │ 34 62 │ 47 65 │ 47 59 │ 32 69 │ 44 51 │ 31 Średni czas światła zielonego: 37.875 Średni czas światła czerwonego: 59.625 Stosunek czasu swiatła zielonego do czerwonego: 0.6352201257861635
Okazuje się, że nasz wynik jest dość zbliżony do rzeczywistości, choć żeby mieć więcej pewności, należałoby wykonać znacznie więcej pomiarów ze względu na nieregularny czas świateł.
Załóżmy, że na autostradzie dwupasmowej zostaje utworzone zwężenie drogi do jednego pasa ruchu (np. z powodu robót drogowych). Skutkiem tego przepustowość się zmniejsza, a ponieważ natężenie na autostradzie jest z reguły duże i najprawdopodobniej większe od przepustowości, to zaczyna tworzyć się korek. Rozpatrzmy taką sytuację na odcinku autostrady A4 od Kątów Wrocławskich do Pietrzykowic (10 106 km). Na polskiej stronie rządowej można znaleźć dane z Generalnego Pomiaru Ruchu 2020/2021 dotyczące średniego dobowego ruchu rocznego na tym odcinku [3].
vehicles1 = DataTable(Liczba_Pojazdów_Na_Dobę = [62776], Motocykle = [60], Samochody_Osobowe = [38901], Samochody_Dostawcze = [7874])
vehicles2 = DataTable(Samochody_Ciężarowe_Bez_Przyczepy = [1219], Samochody_Ciężarowe_Z_Przyczepą = [14571], Autobusy = [151])
show(vehicles1)
println("\n")
show(vehicles2)
1x4 DataTable Liczba_Pojazdów_Na_Dobę │ Motocykle │ Samochody_Osobowe │ Samochody_Dostawcze ────────────────────────┼───────────┼───────────────────┼──────────────────── 62776 │ 60 │ 38901 │ 7874 1x3 DataTable Samochody_Ciężarowe_Bez_Przyczepy │ Samochody_Ciężarowe_Z_Przyczepą │ Autobusy ──────────────────────────────────┼─────────────────────────────────┼───────── 1219 │ 14571 │ 151
Z poprzednich pomiarów do obliczeń wykorzystamy średnią długość samochodu $r = 3,58 \mathrm{m}$ oraz prędkość fali drogowej $ v_f = 6,3 \frac{\mathrm{m}}{\mathrm{s}} $, ponieważ wartości te są dość uniwersalne. Na autostradzie w znacznej ilości pojawiają się także inne pojazdy, więc będziemy liczyć je w następujący sposób:
Jeżeli teraz policzymy średnią długość pojazdu, dostaniemy $r_{śr} = 4,68\mathrm{m}$.
r_śr = (60 * r + 38901 * r + 7874 * 1.5 * r + 1219 * 1.5 * r + 14571 * 2 * r + 151 * 2 * r) / 62776
print("Srednia długość pojazdu: rₛᵣ = ", r_śr)
Srednia długość pojazdu: rₛᵣ = 4.681578091446679
Przyjmijmy $v_{śr} = 120 \frac{\mathrm{km}}{\mathrm{h}} = 33,33\frac{\mathrm{m}}{\mathrm{s}} $.
$t = 24 \mathrm{h} = 86400\mathrm{s}$.
Policzmy teraz gęstość:
G = 62776 * r_śr / (33.33 * 86400)
print("Gęstość: g = ", G)
Gęstość: g = 0.10205560357030727
Zakładamy, że gestość jest stała i równa $g$. Podobnie jest w przypadku przepustowości. Przyjmijmy, że przez zwężenie przejeżdża jeden pojazd na 2 sekundy, czyli:
P = 1 * r_śr / 2
print("Przepustowość: p = ", P)
Przepustowość: p = 2.3407890457233393
Mając już wszystkie potrzebne dane, podstawiamy do pierwszego modelu.
dL_1 = (33.33 * G - P) / (1 - G)
print("L' = ", dL_1)
L' = 1.1812805174713894
$L'$ jest stałe, a więc $L(t)$ będzie funkcją liniową: $L(t) \approx 1,181t$.
Podobnie będzie dla drugiego modelu:
dL_2 = (33.33 * G - P) / (1 - G - P/v_f)
print("L' = ", dL_2)
L' = 2.0151090694329463
Zatem $ L(t) \approx 2,015t $.
Zobaczmy, jak wyglądają wykresy obydwu rozwiazań dla jednej godziny.
plot(t -> t, t -> dL_1 * t, 1:3600, xlabel="t", ylabel="L", label="1. model")
display(plot!(t -> t, t -> dL_2 * t, 1:3600, label="2. model"))
println("Długość korka po godzinie:\n- 1. model: ", dL_1*3600, "\n- 2. model: ", dL_2*3600)
Długość korka po godzinie: - 1. model: 4252.609862897002 - 2. model: 7254.392649958607
Widzimy, że korek w pierwszym modelu osiągnął długość ponad 4km, a w drugim ponad 7km, czyli znacznie więcej. Długość korka według drugiego modelu jest 1,69 razy większa od tej otrzymanej z pierwszego modelu. Zauważamy zatem, jak duży wpływ na długość korka na autostradzie ma uwzględnione w drugim modelu opóźnienie.
Wyprowadziliśmy dwa równania, które następnie, po dopasowaniu do prawdziwej sytuacji, rozwiązaliśmy za pomocą metod numerycznych. Wyniki okazały się być zbliżone do rzeczywistości, co potwierdza, że przedstawione modele mogą mieć zastosowanie w modelowaniu korków drogowych. W szczególności model drugi, który uwzględnia opóźnienie w ruszaniu kolejnych samochodów. Z przeprowadzonych badań dowiedzieliśmy się także, jak ogromny wpływ na tworzenie się korka ma choćby niewielka zmiana czasów świateł zielonego i czerwonego na skrzyżowaniu z sygnalizacją.